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ABSTRACT 

These lectures describe the implementation of optimiza- 
tion techniques based on control theory for airfoil and 
wing design. In previous studies [10, 11] it was shown 
that control theory could be used to devise an effective 
optimization procedure for two-dimensional profiles in 
which the shape is determined by a conformal transfor- 
mation from a unit circle, and the control is the mapping 
function. Recently the method has been implemented in 
an alternative formulation which does not depend on con- 
formal mapping, so that it can more easily be extended 
to treat general configurations [16]. The method has also 
been extended to treat the Euler equations, and results 
are presented for both two and three dimensional cases, 
including the optimization of a swept wing. 

1 FORMULATION OF THE DESIGN PROBLEM 

AS A CONTROL PROBLEM 

Ultimately, the designer seeks to optimize the geometric 
shape of a configuration taking into account the trade-offs 
between aerodynamic performance, structure weight, and 
the requirement for internal volume to contain fuel and 
payload. The subtlety and complexity of fluid flow is 
such that it is unlikely that repeated trials in an interactive 
analysis and design procedure can lead to a truly opti- 
mum design. Progress toward automatic design has been 
restricted by the extreme computing costs that might be 
incurred from brute force numerical optimization. How- 
ever, useful design methods have been devised for vari- 
ous simplified cases, such as two-dimensional airfoils in 
viscous flows [17] and wings in inviscid flows. The com- 
putational costs for these methods result directly from the 
vast number of flow solutions that are required to obtain 
a converged design. 

Alternatively, it has been recognized that the designer 
generally has an idea of the kind of pressure distribu- 
tion that will lead to the desired performance. Thus, it is 
useful to consider the inverse problem of calculating the 
shape that will lead to a given pressure distribution. The 
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method is advantageous, since only one flow solution is 
required to obtain the desired design. Unfortunately, a 
physically realizable shape may not necessarily exist, un- 
less the pressure distribution satisfies certain constraints. 
Thus the problem must be very carefully formulated. 

The problem of designing a two-dimensional profile 
to attain a desired pressure distribution was first studied 
by Lighthill, who solved it for the case of incompressible 
flow with a conformal mapping of the profile to a unit 
circle [13]. The speed over the profile is 

9 =^ 1 , 

where <j> is the potential which is known for incompress- 
ible flow and h is the modulus of the mapping function. 
The surface value of h can be obtained by setting q = qj, 
where q d is the desired speed, and since the mapping func- 
tion is analytic, it is uniquely determined by the value of 
h on the boundary. A solution exists for a given speed q ^ 
at infinity only if 



and there are additional constraints on q if the profile is 
required to be closed. 

The difficulty that the objective may be unattainable 
can be circumvented by regarding the design problem as 
a control problem in which the control is the shape of the 
boundary. A variety of alternative formulations of the 
design problem can then be treated systematically within 
the framework of the mathematical theory for control of 
systems governed by partial differential equations [14]. 
This approach to optimal aerodynamic design was intro- 
duced by Jameson [10, 11], who examined the design 
problem for compressible flow with shock waves, and 
devised adjoint equations to determine the gradient for 
both potential flow and also flows governed by the Euler 
equations. More recently Ta’asan, Kuruvila, and Salas, 
implemented a one shot approach in which the constraint 
represented by the flow equations is only required to be 
satisfied by the final converged solution [20]. Pironneau 
has also studied the use of control theory for optimum 
shape design of systems governed by elliptic equations 
[15]. 
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Suppose that the control is defined by a function 7(() 
of some independent variable £ or in the discrete case a 
vector with components 7 { . Also suppose that the desired 
objective is measured by a cost function I. This may, for 
example, measure the deviation from a desired surface 
pressure distribution, but it can also represent other mea- 
sures of performance such as lift and drag. Thus the 
design problem is recast into a numerical optimization 
procedure. This has the advantage that if the objective, 
say, of a target pressure distribution, is unattainable, it is 
still possible to find a minimum of the cost function. Now 
a variation 67 in the control produces a variation 61 in 
the cost. Following control theory, SI can be expressed 
to first order as an inner product 


requires a number of additional flow calculations equal to 
the number of design variables. Using control theory, the 
governing equations of the flowfield are introduced as a 
constraint in such a way that the final expression for the 
gradient does not require reevaluation of the flowfield. In 
order to achieve this bw must be eliminated from (2). The 
governing equation R expresses the dependence of w and 
7 within the flowfield domain D, 


R (w , 7) — 0, 


Thus 6w is determined from the equation 


SR = 


'dR' 

dw 


bw + 



67 = 0. 


( 3 ) 


si = (g,sr>, 

where the gradient Q is independent of the particular vari- 
ation bT, and can be determined by solving an adjoint 
equation. For a discrete system of equations 


{Q y bT) = Y^ Q ' hTi 


and for an infinitely dimensional system 
{G,6T) = j G(06?d£. 
In either case, if one makes a shape change 


6T = -XQ, (1) 

where A is sufficiently small and positive, then 
61 = - HG,G)<0 
assuring a reduction in I. 

For flow about an airfoil or wing, the aerodynamic 
properties which define the cost function are functions of 
the flow-field variables ( w ) and the physical location of 
the boundary, which may be represented by the function 
T, say. Then 

I = I(w,T) , 

and a change in T results in a change 

< 2 > 

in the cost function. As pointed out by Baysal and Ele- 
shaky [2] each term in (2), except for bw, can be easily 
obtained. and can be obtained directly without 
a flowfield evaluation since they are partial derivatives. 
bT can be determined by either working out the exact 
analytical values from a mapping, or by successive grid 
generation for each design variable, so long as this cost 
is significantly less then the cost of the flow solution. 
Brute force methods evaluate the gradient by making a 
small change in each design variable separately, and then 
recalculate both the grid and flow-field variables. This 


Next, introducing a Lagrange Multiplier -0, we have 

- {£-^[saM 


61 = 


di T 
d7 


dR 1 
37] } 


67 


Choosing to satisfy the adjoint equation 

■Mf. di_ 

dw\ dw 

the first term is eliminated, and we find that 


( 4 ) 


bl = GbT 


( 5 ) 


where 



W 
dT ' 


The advantage is that (S) is independent of bw, with the 
result that the gradient of I with respect to an arbitrary 
number of design variables can be determined without 
the need for additional flow-field evaluations. The main 
cost is in solving the adjoint equation (4). In general, the 
adjoint problem is about as complex as a flow solution. If 
the number of design variables is large, the cost differen- 
tial between one adjoint solution and the large number of 
flowfield evaluations required to determine the gradient 
by brute force becomes compelling. Instead of introduc- 
ing a Lagrange multiplier, V>, one can solve (3) for 6w 



and insert the result in (2). This is the implicit gradient 
approach, which is essentially equivalent to the control 
theory approach, as has been pointed out by Shubin and 
Frank [18, 19]. In any event there is an advantage in 
determining the gradient Q by the solution of the adjoint 
equation. 

After making such a modification, the gradient can be 
recalculated and the process repeated to follow a path of 
steepest descent (1) until a minimum is reached. In order 
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to avoid violating constraints, such as a minimum accept- 
able wing thickness, the gradient may be projected into 
the allowable subspace within which the constraints are 
satisfied. In this way one can devise procedures which 
must necessarily converge at least to a local minimum, 
and which can be accelerated by the use of more so- 
phisticated descent methods such as conjugate gradient or 
quasi-Newton algorithms. There is the possibility of more 
than one local minimum, but in any case the method will 
lead to an improvement over the original design. Fur- 
thermore, unlike the traditional inverse algorithms, any 
measure of performance can be used as the cost function. 

The next section presents the formulation for the case 
of airfoils in transonic flow. The governing equation is 
taken to be the transonic potential flow equation, and the 
profile is generated by conformal mapping from a unit 
circle. Thus the control is taken to be the modulus of 
the mapping function on the boundary. This leads to a 
generalization of Lighthill’s method both to compressible 
flow, and to design for more general criteria. Numeri- 
cal results are presented in Section 3. The mathematical 
development resembles, in certain respects, the method 
of calculating transonic potential flow developed by Bris- 
teau, Pironneau, Glowinski, Periaux, Perrier and Poirier, 
who reformulated the solution of the flow equations as a 
least squares problem in control theory [3]. 


where the density is given by 


P - 




(i 



fptt 


00 


while 



( 8 ) 


Here Moo is the Mach number in the free stream, and the 
units have been chosen so that p and q have a value of 
unity in the far field. 

Suppose that the domain D exterior to the profile C 
in the 2 -plane is conformally mapped on to the domain 
exterior to a unit circle in the a -plane as sketched in 
Figure 1 . Let R and 9 be polar coordinates in the <r-plane, 
and let r be the inverted radial coordinate Also let h 
be the modulus of the derivative of the mapping function 


A = 


dz 

da 


( 9 ) 


Now the potential flow equation becomes 

^(j><f>e) + r^(rp<t> r ) = 0 in D, (10) 

where the density is given by equation (7), and the cir- 
cumferential and radial velocity components are 


2 AIRFOIL DESIGN FOR POTENTIAL FLOW US- 
ING CONFORMAL MAPPING 



(11) 


Consider the case of two-dimensional compressible invis- 
cid flow. In the absence of shock waves, an initially irro- 
tational flow will remain irrotational, and we can assume 
that the velocity vector q is the gradient of a potential <j>. 
In the presence of weak shock waves this remains a fairly 
good approximation. 


D 

c 



la: ^-Plane. 1b: <r-Plane. 

Figure 1 : Conformal Mapping. 

Let p, p, c, and M be the pressure, density, speed-of- 
sound, and Mach number qjc . Then the potential flow 
equation is 

V- = 0, (6) 


while 

q 2 = u 2 + v 2 . (12) 

The condition of flow tangency leads to the Neumann 
boundary condition 

d = -1 ^ = 0 on C. (13) 

h or 

In the far field, the potential is given by an asymptotic 
estimate, leading to a Dirichlet boundary condition at r = 
0[6], 

Suppose that it is desired to achieve a specified veloc- 
ity distribution qj on C. Introduce the cost function 

/ = i £ (q - q d f dO, 

The design problem is now treated as a control problem 
where the control function is the mapping modulus h, 
which is to be chosen to minimiz e / subject to the con- 
straints defined by the flow equations (6-13). 

A modification 6h to the mapping modulus will result 
in variations 6<f>, 6u, 6 v, and 6p to the potential, velocity 
components, and density. The resulting variation in the 
cost will be 

61 = I (q-qd)6qd0, (14) 

Jc 
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where, on C,q = u. Also, 
6u = r 


— 2 6<f) r 

u — , 6v = r — 


6h 

h 


6h 

’~h' 


while according to equation (7) 

dp _ pu dp _ pv 

du c 2 ’ dv c 2 

It follows that 6<j> satisfies 


L6<t> = - 


d_ 

dO 


( 'pM 2 d>e 


) 


~ r Jr (? M2r *’T) 


second term is deleted the method reduces to a variation 
of Lighthill ’s method [13]. 

Equation (19) can be further simplified to represent 
51 purely as a boundary integral because the mapping 
function is fully determined by the value of its modulus 
on the boundary. Set 


log^ = 


where 


^=log 


dz 

da 


= log*, 


where 


L 



puv d \ 

d puv d 1 
dr c 2 30 / 


.(15) 


and 



Then T satisfies Laplace’s equation 


&F = 0 in D } 


Then, if is any periodic differentiable function which 
vanishes in the far field, 

/ 4 U><j>dS= [ P M 2 V<t>-Vxl> 6 -^dS, (16) 
Jd r Jd « 


and if there is no stretching in the far field, T 0. 
Also 5 T satisfies the same conditions. Introduce another 
auxiliary function P which satisfies 

A P = pM 2 VV-VV> in D, (20) 


where dS is the area element r dr d6 , and the right hand 
side has been integrated by parts. 

Now we can augment equation ( 1 4) by subtracting the 
constraint (16). The auxiliary function rp then plays the 
role of a Lagrange multiplier. Thus, 


11 - 

- t %L6tl>dS+ [ pM 2 V4>-Vi> 8 -^ dS. 
Jd r Jd h 


Now suppose that satisfies the adjoint equation 


Lip = 0 in D 


(17) 


with the boundary condition 

dtp_ _ l d f q-qd \ r 

dr pd0\ h ) C ' 

Then, integrating by parts, 

/ %L6<i>dS=- / pxl>MdO, 
Jd r Jc 


and 


61 = - /(«- qd)q^d0 

+ [ pM 2 V<t>-Vxl> 6 -^ dS. 
Jd « 


(18) 


(19) 


Here the first term represents the direct effect of the 
change in the metric, while the area integral represents 
a correction for the effect of compressibility. When the 


and 

P = 0 on C. 

Then, the area integral in equation (19) is 

[ A P6TdS= [ 6F^-d0r [ PA6PdS , 
Jd Jc dr J D 

and finally 


61 


= J Q6J=dO , 


where T c is the boundary value of T, and 
Q - ^ -(.q -q<i)q. 


( 21 ) 


This suggests setting 


6T C = -A G 

so that if A is a sufficiently small positive quantity 

61 = - f A G 2 dO < 0. 

Jc 

Arbitrary variations in T cannot, however, be admitted. 
The condition that T -+ 0 in the far field, and also the 
requirement that the profile should be closed, imply con- 
straints which must be satisfied by T on the boundary C. 
Suppose that log (j^) is expanded as a power series 
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where only negative powers are retained, because other- 
wise (—) would become unbounded for large a. The 
condition that T — ► 0 as a — ► oo implies 

C c = 0 . 

Also, the change in z on integration around a circuit is 
f dz 

Az = / — da = 2tt 2 ci , 

7 aor 

so the profile will be closed only if 

ci = 0. 

In order to satisfy these constraints, we can project Q onto 
the admissible subspace for T c by setting 

c 0 = ci = 0. (23) 

Then the projected gradient G is orthogonal to Q — Q, and 
if we take 

= -A Q y 

it follows that to first order 

si = - f xggde = - f x(g + g-g)gdo 

Jc Jc 

= - f xg 2 de< o. 

Jc 

If the flow is subsonic, this procedure should converge 
toward the desired speed distribution since the solution 
will remain smooth, and no unbounded derivatives will 
appear. If, however, the flow is transonic, one must allow 
for the appearance of shock waves in the trial solutions, 
even if qd is smooth. Then q- qd is not differentiable. This 
difficulty can be circumvented by a more sophisticated 
choice of the cost function. Consider the choice 

7 =U( a,2!+a ’ (£)’)*' < 24 > 

where Aj and A 2 are parameters, and the periodic function 
Z(6) satisfies the equation 


the formula for the gradient becomes 

dP 

G = ^--Zq (27) 

instead of equation (21). Smoothing can also be intro- 
duced directly in the descent procedure by choosing S T c 
to satisfy 

6T '~y **•=-*• < 28 > 

where 8 is a smoothing parameter. Then to first order 

The smoothed correction should now be projected onto 
the admissable subspace. 

The final design procedure is thus as follows. Choose 
an initial profile and corresponding mapping function T. 
Then: 

1. Solve the flow equations (6-13) for <t>, tz, u, q , p. 

2. Solve the ordinary differential equation (25) for Z. 

3. Solve the adjoint equation (15 and 17) or jp subject 
to the boundary condition (26). 

4. Solve the auxiliary Poisson equation (20) for P. 

5. Evaluate G by equation (27) 

6. Correa the boundary mapping function T c by 
calculated from equation (28), projected onto the 
admissable subspace defined by (23). 

7. Return to step 1. 

3 NUMERICAL TESTS OF OPTIMAL AIRFOIL 
DESIGN FOR POTENTIAL FLOW USING CON- 
FORMAL MAPPING 


&Z 

A ' 2 ” X2 ~d&* = q ~ 9d ' 


(25) 


Then, 


61 = 


f c (x,z t z + x,§±iiz)* 

-J c z ( x ' s *- x 4 >z ) m =J c 


Z6qd6. 


Thus, Z replaces q - qd in the previous formulas, and if 
one modifies the boundary condition (18) to 


dip 


\_d_ 

p 86 



on C, 


(26) 


The practical realization of the design procedure depends 
on the availability of sufficiently fast and accurate numer- 
ical procedures for the implementation of the essential 
steps, in particular the solution of both the flow and the 
adjoint equations. If the numerical procedures are not 
accurate enough, the resulting errors in the gradient may 
impair or prevent the convergence of the descent proce- 
dure. If the procedures are too slow, the cumulative com- 
puting time may become excessive. In this case, it was 
possible to build the design procedure around the author ’s 
computer program FL036, which solves the transonic 
potential flow equation in conservation form in a domain 
mapped to the unit disk. The solution is obtained by a very 
rapid multigrid alternating direction method. The original 
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scheme is described in Reference [7]. The program has 
been much improved since it was originally developed, 
and well converged solutions of transonic flows on a mesh 
with 128 cells in the circumferential direction and 32 cells 
in the radial direction are typically obtained in 5-20 multi- 
grid cycles. The scheme uses artificial dissipative terms 
to introduce upwind biasing which simulates the rotated 
difference scheme [6], while preserving the conservation 
form. The alternating direction method is a generalization 
of conventional alternating direction methods, in which 
the scalar parameters are replaced by upwind difference 
operators to produce a scheme which remains stable when 
the type changes from elliptic to hyperbolic as the flow 
becomes locally supersonic [7]. The conformal mapping 
is generated by a power series of the form of equation (22) 
with an additional term 

(■-SM 1 -;) 

to allow for a wedge angle t at the trailing edge. The 
coefficients are determined by an iterative process with 
the aid of fast Fourier transforms [6], 

The adjoint equation has a form very similar to the 
flow equation. While it is linear in its dependent variable, 
it also changes type from elliptic in subsonic zones of the 
flow to hyperbolic in supersonic zones of the flow. Thus, 
it was possible to adapt exactly the same algorithm to 
solve both the adjoint and the flow equations, but with re- 
verse biasing of the difference operators in the downwind 
direction in the adjoint equation, corresponding to the re- 
versed direction of the zone of dependence. The Poisson 
equation (20) is solved by the Buneman algorithm. 

An alternative procedure would be to derive the ex- 
act adjoint equation corresponding to the discrete equa- 
tions which approximate the potential flow equation. This 
would produce the exact derivative of the discrete cost 
function with respect to the discrete control, at the ex- 
pense of very complicated formulas and a costly inversion 
procedure. The discrete adjoint equation would then be 
a particular discretization of the differential adjoint equa- 
tion corresponding precisely to the discretization used for 
the flow equation. The efficiency of the present approach, 
which uses separate discretizations of the flow and ad- 
joint equations, depends on the fact that in the limit of 
zero mesh width the discrete adjoint solution converges 
to the true adjoint solution. This allows the use of a rather 
simple discretization of the adjoint equation modeled after 
the discretization of the flow equation. Numerical exper- 
iments confirm that in practice separate discretizations of 
the flow and adjoint equations yields good convergence 
to an optimum solution. 

As an example of the application of the method. Fig- 
ure 3 presents a calculation in which an airfoil was re- 
designed to improve its transonic performance by reduc- 
ing the pressure drag induced by the appearance of a shock 
wave. The drag coefficient was therefore included in the 


cost function so that equation (24) is replaced by 

7= ix( Aiz2+Aj ( d 4y) M+x ’ c - 

where A 3 is a parameter which may be varied to alter the 
trade-off between drag reduction and deviation from the 
desired pressure distribution. Representing the drag as 

D = J c (p - p ” ) te il> ' 

the procedure of Section 2 may be used to determine the 
gradient by solving the adjoint equation with a modified 
boundary condition. A penalty on the desired pressure 
distribution is still needed to avoid a situation in which 
the optimum shape is a flat plate with no lift and no drag. 

It was also desired to preserve the subsonic charac- 
teristics of the airfoil. Therefore two design points were 
specified, Mach 0.20 and Mach 0.720, and in each case 
the lift coefficient was forced to be 0.6. The composite 
cost function was taken to be the sum of the values of the 
cost function at the two design points. The transonic drag 
coefficient was reduced from 0.01 91 toO.OOOl in8design 
cycles. In order to achieve this reduction the airfoil had 
to be modified so that its subsonic pressure distribution 
became more peaky at the leading edge. This is consis- 
tent with the results of experimental research on transonic 
airfoils, in which it has generally been found necessary 
to have a peaky subsonic presure distribution in order to 
delay the onset of the transonic drag rise. It is also impor- 
tant to control the adverse pressure gradient on the rear 
upper surface, which can lead to premature separation of 
the viscous boundary layer. It can be seen that there is no 
steepening of this gradient due to the redesign. 

4 DESIGN FOR POTENTIAL FLOW USING A FI- 
NITE VOLUME DISCRETIZATION SCHEME 

While the use of conformal mapping, as it has been pre- 
sented in sections 2 and 3, leads to an effective design 
method for two dimensional profiles, it is not easy to treat 
more complex configurations because of the difficulty in 
devising appropriate numerical mapping methods. More- 
over, conformal mapping is limited to two dimensional 
transformations. In this section an alternative formula- 
tion using a general coordinate transformation is adopted. 
This is intended to be a precursor to the three dimensional 
problem. 

Consider the case of two-dimensional compressible 
inviscid flow. A general transformation from cartesian 
coordinates x and y to the coordinates f and r\ can be 
represented by the transformation 


ax 

dx 

n 

dv 



L ot 

dt) J 


6 


The potential flow equation can be written in divergence 
form as 


A A 

dx + dy ^ = ° ' nD ' ^ 

where u and v represent the Cartesian velocity compo- 
nents. The coordinate transformations may be defined 


CHU ■ 


dx dx 
6i aa 

L d y °v 


h 

4> 


:} 


= K 1 


i::}- 


(30) 


Also 


/ \ _ If { ** \ = K T { ^ \ 


Then 


jz(pJU)+ 7r (pJV) = 0 in D. (31) 
dt Otj 


where J is the Jacobian 


•'=-**>= I™- 


defined by the flow equations (29-33). The first variation 
of the cost function is 

81 = J 
- 

+ l « -«>${%)%■*■ o5 > 

since on the wall 

_<W_ d±di 

qw ~ ds~ dids ' 

In general we need to find how a modification to the airfoil 
geometry causes a variation 8<j> y as well a variation in the 
grid parameters 6 An, 8 A\ 2 , 8 A 22 , and 6J . The variations 
in U,V and p are 

SU = 6 (An) <t>( + An<5<^ A S (A 12 ) (f>r) + Ai 28 <fir) 

S V = 8 (A 12 ) <t>{ 4- A\ 28 <f>£ -h 8 (A 22 ) <f*r) + A22<5<£tj 


Here, U and V represent the contravariant velocities 



Thus, 

U = An<j>{ + Ai2<f> v (32) 

V = Ai2tt + A2rt n . (33) 


where 


A = (K T K)~ ] = 


An 

A\2 


A \2 

A 22 


Consider first the case in which the cost function is 
defined such as to achieve a target speed distribution: 


I = 



(34) 


where q * is the desired speed distribution and C is the 
airfoil surface. 

The design problem is now treated as a control prob- 
lem where the control function is the airfoil shape, which 
is to be chosen to minimize I subject to the constraints 


8 P = -3 


vUvP 


c 2 [ d£ dr) J 


6<j> 


[ti-An<l>£ + 26Ai2$ n fy +6 Az24%] ■ 


2c 2 

It follows that 6<t> satisfies 


L6<f> - 


where 


L = 


-■^-Q(SJ,6An,6Ai2, 6A22) 

d< 

—• 3 - P(6J , 6A\\,f>A\2, &A 22 ), 

dr) 


d_ 

dt 


pj (An - ^r) 

+pJ (A\2 — — 


6_ 

8i 

) h J 


+ 7T“ 


uv\ a_ i 
tT I ae 


Q pj (A 12 

dr) +pJ (A 22 - tt) J 


(36) 


(37) 


and 


Q(6J,6A n ,6A 12 ,6A22) = pUSJ 

+pJ<t>i - -^f) 

+pj+, (l - U -jf) M» 

+ e ] ** (~ W ) iM ' 

P(6J,6A U ,6A 12 ,SA 2 2) = pV6J 
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V+n\ c A 

” SA 22 


+pJ <Pr) ( 1 — 


+p j 4( ( 1 - 


) < 5-^12 
+pj<j>( (— 5 ^) 6 A U- 


If tp is any periodic function vanishing in the far field* 
equation (36) can be multiplied by ip and integrated over 
the domain. After integrating the right hand side by parts 
we arrive at 

J.* 1 ***" 

+ [ {rPpJ [ 6 A l 2 <j> ( + } d(. (38) 

Jc 

Now subtracting (38) from (35), 
d (q - qa) 


61 


-L 


6 <t>d£ 


i c dt 
f 1 i ( ds\ 

+ J c 2 «- t * S {w« 

^ [ f .d<i> x {dt\ds, e 

— J ipL 6 <pd£dri 

+ l% Q+ !± Pd(dv 


V 


+ 


f {iPpJ [f>Ai2<p{ + 6A22<Ptj] } d£. 

Jc 


(39) 


Then setting up the adjoint system we have 
Lip = 0 in D, 
with the boundary condition 

pj + A-n^) = — (q - qa) ( 40 ) 

After applying the second form of Green ’s theorem to (39) 
we get 

f il>L 6 <t>dS = f {ippJ [ 6 A\ 2 <i>t J r&A 22 <i > n \ } d£ 
Jd Jc 

+ J { 6 <f>pJ + 6 A 22 ipri ] } d£. 

Finally the variation can be defined as 

■L*-» *( 2 )S« 

<4i) 


+ 

+ 


No general analytic grid transformation is generally 
available for the finite volume formulation. Furthermore, 
the variation with respect to the grid quantities is now 
spread into 6 An, 6 A i2 , 6 A 2 2, and 6 J instead of just the 
modulus of the transformation as was the case for confor- 
mal mapping. Therefore, to construct <57, an independent 
basis space of perturbation functions fc< , i = 1 , 2 , . . . , n 
(n = number of design variables) is chosen that allows 
for the needed freedom of the design space. Thus, the 
shape T now becomes !F(bi), where the functions 6, now 
represent the control. The variations 6 A U , 6 A \ 2 , 6 A 22 , 
and 6 J are obtained by a direct finite difference proce- 
dure with respect to each design variable 6<. Once 61 is 
obtained, any optimization procedure can be employed to 
minimize the cost with respect to the given basis 6, . 

If the flow is subsonic, this procedure should converge 
toward the desired speed distribution since the solution 
will remain smooth, and no unbounded derivatives will 
appear. If, however, the flow is transonic, one must allow 
for the appearance of shock waves in the trial solutions, 
even if qa is smooth. In such instances q — qa is not dif- 
ferentiable. As in section 2, the cost function is redefined 
as 


-u( 


X\Z 2 + x 2 



where Ai and A 2 are parameters, and the periodic function 
£(£) satisfies the equation 


d?Z 

X\Z - A 2 -^J = q - q d . 

Then, 

61 = iy zsz+ ^i sz ) 


(42) 


di 


= J c Z(x l 6Z-X 2 J^6Z S ) d£ = J Z6qdt. 

Thus, Z replaces q — qa in the previous formula and one 
modifies the boundary condition (40) to 

pJ (AiaV’f + A 22 <t>ti) = — (Z) on C. (43) 

For the case where the cost function is drag, (34) is 
replaced by, 

tm L' %*■ <44 > 

The first variation of the cost function is now, 

w - jy^ym- 

■ *«)« 

*M 8 


di. 


(45) 
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Thus, (41) becomes 


61 = - 


+ 

+ 


Jc P6 (W d< 


(46) 


where the boundary condition on ip y (40) or (42), is re- 
placed with 


pj + ^ 22 <Prf) 



(47) 


or 


XiZ 




(48) 


The entire procedure can be summarized for the cost func- 
tion based on target speed distribution as follows: 


1. Solve the flow equations (29-33) for </>, u y v y q, p y 
U y and V. 

2. Smooth the cost function if necessary by (42). 

3. Solve the adjoint equation (37 and 39) for ip subject 
to the boundary condition (40) or (43). 

4. For each t independently perturb the design vari- 
ables, &,•, and calculate the necessary metric vari- 
ations ( 6 A Uy 6 A\ 2 » 6 A 22 , <57, and *(&)) byre- 
calculating the perturbed grid with automatic grid 
generation. 

5. Directly evaluate 61 by equation (41). 

6. Project 61 into a feasible direction subject to any 
constraints to obtain 6L 

7. Feed 6 1 as the gradient with respect to 6, to a quasi- 
Newton optimization procedure. 

8. Calculate the search direction with a quasi-Newton 
algorithm and perform a line search. 

9. Return to 1 if the process has not converged. 


In practice the method resembles that used by Hicks et 
al [17] with the control theory replacing the brute force, 
finite difference based, gradient calculation. The current 
formulation has an advantage by requiring computational 
work proportional to 2 + m flow solver evaluations (m 
being the number of calculations required per line search) 
per design cycle as opposed to 1 + m + n. Thus, un- 
like conventional design optimization programs, the cur- 
rent method’s computational cost does not hinge upon the 
number of design variables provided the grid regeneration 
is fast and automatic. The method also has the advantage 
of being quite general, allowing arbitrary choices for both 
the design variables and the optimization technique. 


5 NUMERICAL IMPLEMENTATION OF THE 
GENERALIZED POTENTIAL FLOW DESIGN 
METHOD 

The practical implementation of the generalized poten- 
tial flow design method, as with the conformal potential 
method, relies heavily upon fast accurate solvers for both 
the state (<p) and co-state (ip) fields. Further, to improve 
the speed and realizability of the methods, a robust choice 
of the optimization algorithm must be made. Finally, ap- 
propriate design variables must be chosen which allow 
sufficient freedom in realizable designs. In this work, the 
author’s FL042 full potential computer program and the 
QNMDIF (by Gill, Murray and Wright [4]) quasi-Newton 
optimization algorithm are employed. 

In FL042 the flow solution is obtained by a rapid 
multigrid alternating direction method [7]. The scheme 
uses artificial dissipative terms to introduce upwind bi- 
asing which simulates the rotated difference scheme [6] 
while preserving the conservation form. The alternating 
direction method is a generalization of conventional alter- 
nating direction methods in which the scalar parameters 
are replaced by upwind difference operators to produce 
a scheme which remains stable as the equations change 
type from elliptic to hyperbolic in accordance with the 
flow becoming locally supersonic [7], 

QNMDIF is an unconstrained quasi-Newton optimiza- 
tion algorithm that calculates updates to a Cholesky fac- 
tored Hessian matrix by the BFGS (Broyden-Fletcher- 
Goldfarb-Shanno) rank-two procedure. Hence, informa- 
tion about the curvature of the design space feeds in 
through the successive gradient calculations. 

Since the primary computational costs arise from not 
only the flow solution algorithm but also the adjoint solu- 
tion algorithm, both need to be computationally efficient. 
The adjoint equation has a form very similar to the flow 
equation. While it is linear in its dependent variable, 
it also changes type from elliptic (in subsonic zones of 
the flow) to hyperbolic (in supersonic zones of the flow). 
Thus, it was possible to adapt exactly the same algorithm 
to solve both the adjoint and the flow equations, but with 
reverse biasing of the difference operators in the down- 
wind direction for the adjoint equation, corresponding 
to its reversed direction of the zone of dependence. A 
multigrid method is used to accelerate the convergence 
of a generalized alterating direction scheme in a manner 
similar to the flow solver. 

Design variables are chosen with the following form, 
suggested by Hicks and Henne [5]: 

b(x) = sin ( 7rx kftoOi) J 

b(x) = x tl (1 — x)e“ <2 *, 

where t \ and 1 2 control the center and thickness of the per- 
turbation and x is the normalized chord length. When dis- 
tributed over the entire chord on both upper and lower sur- 
faces these analytic perturbation functions admit a large 
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possible design space. They have the advantage of be- 
ing space based functions, as opposed to frequency based 
functions, and thus they allow for local control of the de- 
sign. They can be chosen such that symmetry, thickness, 
or volume can be explicitly constrained. Further, particu- 
lar choices of these variables will concentrate the design 
effort in regions where refinement is needed, while leav- 
ing the rest of the airfoil section virtually undisturbed. The 
disadvantage of these functions is that they do not form a 
complete basis space, nor are they orthogonal. Thus, they 
do not guarantee that a solution, for example, of the in- 
verse problem for a realizable target pressure distribution 
will necessarily be attained. Here they are employed due 
to their ease of use and ability to produce a wide variation 
of shapes with a limited number of design variables. 

The generalized potential flow design algorithm based 
on the finite volume scheme has been applied to a vari- 
ety of test cases, which are described in the following 
paragraphs. These include both non-lifting cases, where 
a symmetric taiget pressure distribution is specified and 
the optimization is started from an arbitrary symmetric 
initial guess, and lifting cases where the target pressure 
distribution is specified, and finally cases which verify the 
capability of the method to find profiles with minimum 
drag. 

The first non-lifting example shown in Figure 4, illus- 
trates that for subsonic flow, M oo = 0.2 and a = 0° , a 
given airfoil shape, in this case a NACA 64012, can be 
recovered by starting from an arbitrary shape and speci- 
fying the target pressure distribution. A close look at the 
final solution shows that a small discrepancy is evident at 
the trailing edge. This may be associated with the lack of 
completeness of our basis space. In the next example, see 
Figure 5, the design takes place at = 0.8, a = 0°, 
where the initial NACA 0012 airfoil is driven towards the 
subsonic pressure distribution of the NACA 64021. In 
this case the target pressure distribution exceeds Cp* for 
Moo = 0.8. Therefore, the pressure distribution repre- 
sents shock free transonic flow. Since, in general, such a 
pressure distribution may not be realizable, the program 
approaches the taiget with the nearest feasible pressure 
distribution. An examination of Figure 5 demonstrates 
that a very weak shock in the designed pressure distribu- 
tion replaces the smooth transition to subsonic flow seen 
in the target distribution. In the final example non-lifting 
case of Figure 6, an arbitrary pressure distribution which 
does contain a shock wave and is realizable, is used as 
the target. Here the computer program was able to obtain 
the corresponding airfoil geometry along with the correct 
shock wave location with a high degree of accuracy, as 
can be seen both in the pressure distribution and in the 
airfoils. 

The second group of test cases address the problem of 
attaining a desired pressure distribution for lifting airfoils. 
The most convenient method of obtaining such solutions 
with the present design method is to determine the lift co- 
efficient associated with the target pressure distribution, 


and match this lift with the initial airfoil. The design pro- 
gresses with the flow solver and the adjoint system being 
driven by constant circulation instead of fixed angle of 
attack. The first example using this technique, shown in 
Figure 7, drives the NACA 0012 airfoil toward the tar- 
get pressure distribution for the NACA 64A410 airfoil at 
Moo - 0.735, a = 0°,andCj = 0.75. This case requires 
a shift in the shock location and a significant change in 
the profile shape such that the target pressure distribution 
is obtained. The final solution almost exactly recovers 
the pressure distribution and the airfoil shape. In the next 
example, Figure 8, the NACA 0012 airfoil is again used 
as the starting condition to obtain the pressure distribution 
of the GAW72 airfoil operating at Moo = 0.7, a = -2° , 
and Ci = 0.57. This case is difficult since the target air- 
foil has a cusped trailing edge while the initial airfoil has a 
finite trailing edge. As was seen in some of the non-lifting 
cases, there are small discrepancies evident near the trail- 
ing edge that may be due to the incomplete basis of the 
chosen design variables. The difference in the profiles 
between the final design and actual GAW72 is partly due 
to the fact that the GAW72 coordinates place the trailing 
edge at a non-zero y ordinate while the NACA 00 1 2 places 
the trailing edge at y = 0. Also, the redesigned airfoil is 
subject to an arbitrary rotation since the angle of attack is 
free during optimization. The last test case in which the 
design program is run in inverse mode involves driving 
the NACA 0012 airfoil at M w = 0.75 to obtain the target 
pressure distribution of the R AE airfoil at the same Mach 
number, a = 1.0°, and C\ = 0.80. Due to the steep 
favorable pressure gradient at the leading edge upper sur- 
face and the strong shock exhibited (see Figure 9) by the 
RAE airfoil at these conditions this case represents quite 
a difficult test for the program. The method recovers the 
taiget pressure distribution almost exactly. A comparison 
of the profiles reveals that the the designed airfoil has no 
observable differences when overlayed with the original 
airfoil. 

The last group of results introduces drag as the cost 
function. Again the design process is carried out in the 
fixed lift mode. In Figure 10, the first drag minimization 
example, a NACA 0012 is again used as a starting airfoil. 
The design takes place at Moo = 0.75 and C\ - 0.50 
where a strong shock causes considerable wave drag in 
the initial airfoil. To make the problem interesting, the 
optimization is carried out such that symmetry of the de- 
sign is preserved. The final design is a symmetric airfoil 
with an increased maximum thickness that operates at the 
same lift coefficient, but has a reduction in drag from 
Cd — 0.0127 to Cd = 0.0016. In the final test case (see 
Figure 9) the camber distribution is optimized instead of 
thickness distribution. The design starts from a NACA 
64A410 airfoil operating at Moo - 0.75, and C\ - 0.60 
which displays 42 counts of drag according to the potential 
flow calculation. By allowing only changes to the camber 
distribution, a final airfoil is produced which maintains 
C\ — 0.60 but does so with only 4 counts of drag. 
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6 DESIGN OF AIRFOILS USING THE EULER 
EQUATIONS 

This section extends the application of control theory for 
aerodynamic shape optimization to the Euler equations for 
two dimensional flow. Consider the case of compressible 
flow over an airfoil. In the absence of separation and other 
strong viscous effects, the flow is well approximated by 
the Euler equations. In contrast to the previous implemen- 
tations which relied on the isentropic potential equation, 
here strong inviscid shocks are modeled correctly with 
entropy production. Consider the flow in a domain D. 
The profile defines the inner boundary C, while the outer 
boundary B is assumed to be distant from the profile. Let 
p y p y u y v y E and H denote the pressure, density, Cartesian 
velocity components, total energy and total enthalpy. For 
a perfect gas 


with 


p = (ry- \)p$yE- 1 (u 2 + t; 2 ) j , 


(49) 


and 


pH = pE + p, (50) 

where 7 is the ratio of the specific heats. The Euler 
equations may then be written as 

^ + S + |=° ” D ' < 51) 

where x and y are Cartesian coordinates, t is the time 
coordinate and 


W : 



/ = < 


pu 

( pv 

pu 2 + p 

k „ _ J pvu 

puv 

' 9 I pv 2 + p 

„ P UH J 

{ pvH t 


(52) 


Consider a coordinate transformations to computa- 
tional coordinates £ , rj with the transformation matrix 


K = 


dx 

dt 


dx 

dr) 

SJL , 
dr} J 


and the Jacobian 

dx dy dx dy 

d£ dy dr} 

Introduce contravariant velocity components 


J = 


KK'CH 


dy dx_ 

dr) drj 


dz 


dx 

Sf 


The Euler equations can be written as 


dW dF 8G „ . _ 

— = 0 mD, 
dt d(, dr} 


C) 


(53) 


W = J < 


P 

pu 

pv 

pE 





f ' 

pU 


pV 

V 

II 

P Uu + %P 

C5 

II 

pVu+^p > 


pUv + | ip 

pVv+ |jp 


pUH 

< j 


pVH 


( 54 ) 


Assume now that the computational coordinate system 
conforms to the airfoil section in such a way that the 
surface C is represented by rj — 0. Then the flow is 
determined as the steady state solution of the equation 
(54) subject to the flow tangency condition 


F = 0 onC, 


(55) 


At the far field boundary B y conditions are specified for 
incoming waves, while outgoing waves are determined 
by the solution. 

Consider the case of the inverse problem where the 
cost function may be defined as 

/ = j/ c 0 >-w)^=i/ c (p - Pi )2 (!)<«. 

where pd is the desired pressure. The design problem 
is now treated as a control problem where the control 
function is the airfoil shape, which is to be chosen to 
minimize I subject to the constraints defined by the flow 
equations (53-55). A variation in the shape will cause a 
variation 6p in the pressure in addition to a variation in 
the geometry and consequently the variation in the cost 
function becomes 


61 


= J c (p ~ Pd)6p {^) <Ui 
4 / c (p_Pd) ^(i)^- 


(56) 


Since p depends on w through the equation of state 
(51-52), the variation 6p can be determined from the 
variation 6w. Define the Jacobian matrices 


* = * = !:• c t = Y. JK » A r < 57 > 


dw' 


Then the equation for 6w in the steady state becomes 


£(«*•)+ £<«©=<>, 


(58) 
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where 


6F = C^ + i(j^)j + s(j g), 
6G = C 2 6w + 6[J^f + s(j 


£)'• 


Now, multiplying by a vector co-state variable i> and 
integrating over the domain 

36G' 


L 


*‘iw + '-w ]d(d,,=o - 


and if -0 is differentiable this may be integrated by parts 
to give 


f ( d^ T 

Jd \ dt 


• T dxb T \ 

6F + lk 6G ) ^ = 


f (ni il> T 6 F + n 2 rp T 6 G) d£ 

Jb 

+ / (n ^ t 6 F + n 2 ^ T 6 G) d £, 

Jc 

where n t * are the components of a unit vector normal 
to the boundary. No boundary integrals appear in the 
77 direction because the mesh is assumed to be of 0 - 
type, with the result that the solution is periodic in the 
£ coordinate thereby canceling the rj boundary integrals. 
Thus the variation in the cost function may now be written 

+ 5 X»-"*(aO« 

+ J D (w eF+ %- 6G ) d(d ’' 

— [ (n\\l) T 6F + n2ip T 6G) d£ 

Jb 

— f {n\\l> T 6 F 4- ri 2 ip T 6 G) d£. 

Jc 

On the profile n x = 0 and n 2 = — 1. It follows from 
equation (55) that 


6G= J < 


0 

§§*P 

I 




0 


( J &) 


0 


(59) 


Suppose now that rj> is the steady state solution of the 
adjoint equation 


c T —-C T —-0 in D 

«<r °2 - 0 mlJ . 


6( (60) 

At the outer boundary incoming characteristics for ^ cor- 
respond to outgoing characteristics for 6w. Consequently, 
one can choose boundary conditions for ^ such that 

iiixl> T Ci6w = 0. 


Then if the coordinate transformation is such that <5 ( JK~ l ) 
is negligible in the far field, the only remaining boundary 
term is 

/ ^ T 6Gdt. 

Jc 

Thus by letting ^ satisfy the boundary condition, 

j(v- 2 g + ^g)=- ( P-Pa ) | onC, (61) 
we find finally that 

= \ J c (P~Pd) 2 6 d£ 

fAA j £) + A j %)}‘ 


61 


+ 


pdC (62) 


If the flow is subsonic, this procedure should converge 
toward the desired pressure distribution since the solution 
will remain smooth, and no unbounded derivatives will 
appear. If, however, the flow is transonic, one must allow 
for the appearance of shock waves in the trial solutions, 
even if pd is smooth. In such instances p - pd is not dif- 
ferentiable. As in the case of potential flow, this difficulty 
can be circumvented by a more sophisticated choice of 
the cost function. Consider the choice 

where Ai and A 2 are parameters, and the periodic function 
Z(£) satisfies the equation 


. 7 . <PZ 

Aj2- A 2 -^j =P~Pd- 


(63) 


Then, 

61 


= L( 

■ SA 
-L 


\\Z 62 + A 2 



\ ds 

) 

P z ) 

p 


dt 


Z6q^di. 

c 


Thus, Z replaces p-pd with a corresponding modification 
to the boundary condition for the adjoint equations. 

A convenient way to treat an airfoil is to use a confor- 
mal mapping of the profile in the z plane to a near circle 
in the a plane, followed by shearing of the radial coor- 
dinate to make the system boundary conforming. Polar 
coordinates are introduced in the mapped plane a. When 
mapped back to the physical plane this gives a smooth. 
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nearly orthogonal grid. This procedure is intermediate 
between the use of a full conformal mapping as in sec- 
tion 2 and an arbitrary numerically generated grid as in 
section 4. We can now specialize our generalized design 
procedure to treat this grid system. Define the first con- 
formal mapping from z to a by letting the derivative of 
the mapping function be 

da 

Now using polar coordinates r, and 0 in the a plane, the 
first transformation matrix is 


XB Xr 

ye y r 


vs c 
— cs s 


and we can define contravariant velocities 


s —c 
c s 


]{:} 


s = sin (/? — 6) , c = cos (/? — 6) . 

The Euler equations can now be represented in the a plane 
as 

^ + ^ + ^ =0 (64) 


P U 

pUu ■+ sp 
pU v — cp 
pUH 


pV 

pVu + cp 
pVv + sp 
pVH 


Now let the final computational coordinates be defined by 
a radial shearing transformation 

6 = t, r = T) + S(t) 

and the transformation matrix 

rn si f i o i 


Q6 

89 


dr) 

dr 

dr 

ae 

d V 


dS 1 
d( 1 


, det(K 2 )=l. 


Now we can identify the complete transformation matrix 
as 

* =*,*,=» i ’■ t+,s<c 


— rc -f S^s s 


while the fluxes are 

hF = h ( sf — eg) 

= Vr,f - Zr)9 

and 

h [{ri + S)G-S ( F] = h[(n + S)c-S(s]f 
+ h [(77 + S)s + <%c] g 
= x ( g-y ( f, 

Thus the Euler equations assume the form 
((*7 + S)h 2 W) 

+ ^(hF) + A (h + S) G - hS^F) = 0, 

while the surface tangency condition on the velocity be- 
comes 

X{V — u = h [(77 + S) V — U] =0. 

Now we take ) as the control. It is also convenient to 
represent the inverse problem by the cost function 

I = ^J c (P~Pd) 2 dd = 1 = l^(p-Pd) 2 ^- 

This eliminates terms in 6 from the gradient. The 

variations in the fluxes become 

6(hF) = C x 6w 

6 [h (rj + S) G - hS{F] = C 2 6 w + hSSG - h 6 S e F 

where C\ and C 2 are the Jacobian matrices defined in 
equation (57). Choosing rp to satisfy the adjoint equation 
(60) with the boundary condition 

X{fo - yrf 2 = h [(t) + S)s + S{c] i/>n=p-pd 

the variation in the cost reduces to 


/ ( P-Pd)6pd£ 

Jc 

[ rl> T (6ShG-6S ( hF)dS 
Jc 

[ ( 6ShG + 6S( hF) d£dT), 

Jd 


where F and G are the fluxes defined in equation (65), 
and F and G are F and G with the pressure terms deleted. 
Define 

P = ip T hF + j ip T -^(hF)drj 
Q = rj> T hG + Jxp T ^-(hG)dg. 
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Then 


61 = J (Q6S-P6Ss)dt 

= [ gssdt, 

Jc 


where the gradient is 


section 2. Thus the change in the shape function S (0 is 
defined by solving 

ls -rPk s= - xe ' 

where (3 is a smoothing parameter. Then, to first order, 
the variation in the cost is 


g = Q+ ik- (66) 

Hie entire procedure can be summarized as follows. 

1. Solve the flow equations (51-55) for p, u, v, p, E, 
H, U, and V. 

2. Smooth the cost function if necessary by (63). 


61 = 


< 


Q6Sdi 


I 

-Jj[ 6s - SS ^M 6S . 




0. 


3. Solve the adjoint equations (60) for ip subject to the 
boundary condition (61). 

4. Calculate P and Q from the variation in the control 

5 ( 0 - 

5. Evaluate Q by equation (66) 

6. Project Q into a feasible direction subject to any 
constraints to obtain C. 

7. Correa the mapping in the direction of steepest 
decent 

6S(Z) = -\G. 


The option to minimize the pressure drag coefficient 


I =C d = 


IP oo?; 


Jc 


dy.. 

p at d t' 

c 


where c is the chord length, has also been included. Tb 
prevent the procedure from trying to reduce drag by re- 
ducing the profile to a non-lifting flat plate a target pres- 
sure distribution is retained in the cost function, which 
becomes 


/= J c (p-p d ) 2 dZ + a 2 C d 


or by using C as the gradient in a quasi-Newton or 
conjugate gradient search method. 

8. Return to 1. 

7 IMPLEMENTATION OF THE EULER BASED 
DESIGN METHOD 

The practical implementation of design method relies 
heavily upon fast accurate solvers for both the state (tu) 
and co-state (ip) fields. Further, to improve the speed 
and realizability of the method, a robust choice of the 
optimization algorithm must be made. In this work, the 
author’s FL082 full computer program has been used to 
solve the Euler equations. This program uses a multi 
stage time stepping scheme with multi grid acceleration 
to obtain very rapid steady state solutions, typically in 25 
steps [8, 9]. The adjoint equations are solved by a similar 
method, in which the flux calculations for the Euler equa- 
tions are replaced by the corresponding formulas for the 
adjoint equation. 

In the initial tests a simple gradient procedure has 
beat used as the optimization process. To preserve the 
smoothness of the profile the gradient is smoothed at each 
step in a similar manner to that used in the method of 


where £2i and £2 2 are weighting parameters. Also the 
calculations are performed at a fixed lift coefficient corre- 
sponding to that of the target pressure distribution, while 
the angle of attack is allowed to vary as needed. Three 
test cases are presented for the design algorithm. The 
first two address the problem of attaining a desired pres- 
sure distribution. The first example using this technique, 
shown in Figure 12, drives the Korn airfoil toward the 
taiga pressure distribution for the NACA 64A410 air- 
foil at Moo — 0.75, a = 0°, and Q = 0.7. This case 
requires a shift in the shock location and a significant 
change in the profile shape such that the targa pressure 
distribution is obtained. The final solution almost exactly 
recovers the pressure distribution and the airfoil shape. In 
the next example. Figure 13, the Kom airfoil operating 
at Moo = 0.78 is used as the starting condition to ob- 
tain the pressure distribution of the same airfoil operating 
at Moo — 0.75, a = 0°, and Ci = 0.64. This case is 
difficult since the targa pressure distribution may not be 
realizable from a physical profile. Note that while the 
achieved pressure distribution is very close to the taiga 
pressure distribution, the drag of 77 counts is much larger 
then the zero drag experienced by the Kom airfoil at its 
design point The third test case introduces drag as the 
cost function. Again the design process is carried out in 
the fixed lift mode. In Figure 14, a NACA 64A410 is 
again used as a starting airfoil. The design takes place at 
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Moo = 0.75 and C\ — 0.68 where a strong shock causes 
considerable wave drag in the initial airfoil. To preserve 
a reasonable lifting airfoil shape the cost function is con- 
structed as a blend of preserving the original pressure 
distribution and reducing the drag. The final design has a 
reduction in drag from Cd — 0.0144 to Cd — 0.0018. 


Assume now that the new computational coordinate sys- 
tem conforms to the wing in such a way that the wing 
surface Bw is represented by £2 = 0. Then the flow is 
determined as the steady state solution of equation (71) 
subject to the flow tangency condition 

— 0 on By/ . (73) 


8 THREE DIMENSIONAL DESIGN USING THE 
EULER EQUATIONS 


In order to illustrate further the application of control 
theory to aerodynamic design problems, sections 8 and 
9 treat the case of three-dimensional wing design, again 
using the inviscid Euler equations as the mathematical 
model for compressible flow. In this case it proves con- 
venient to denote the Cartesian coordinates and velocity 
components by xi, X2, X3 and iq, 1*2, «3, and to use the 
convention that summation over i = 1 to 3 is implied by a 
repeated index i. The three-dimensional Euler equations 
may be written as 


dtv dfj 
dt dxi 


0 in D, 


where 


p 


pui 

pu 1 


puiui + p6n 

i pU2 

► , /< = < 

pUiU 2 + pSi2 > 

PU3 


pUiU 3 + p6i3 

PE J 


„ pUiH 


(67) 


( 68 ) 


and 6ij is the Kronecker delta function. Also, 

P = (7 ~ \)pIe~ 1 (uf)| , (69) 


and 

pH = pE + p (70) 


where 7 is the ratio of the specific heats. Consider a 
transformation to coordinates £1 , £2, £3 where 


\dxA 


' 90 ' 

'IS? 
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■«*> 

£ 

J = det (K) , K t -/ = 

m dxj _ 


Introduce contravariant velocity components as 



The Euler equations can now be written as 

dW dFi A _ 
— - 4- = 0 m Z>, 

dt Ob 


(71) 


with 


r > 

P 


pUi 

pui 


pUM + ffcp 

< pu 2 

> , Fi = J < 

P UiU 2 + ||jp 

pti 3 


pUiti 3 + 

pE 

\ * 




(72) 


At the far field boundary Bp , conditions are specified for 
incoming waves, as in the two-dimensional case, while 
outgoing waves are determined by the solution. 

Suppose now that it is desired to control the surface 
pressure by varying the wing shape. It is convenient 
to retain a fixed computational domain. Variations in 
the shape then result in corresponding variations in the 
mapping derivatives defined by H . Introduce the cost 
function 


I= \JJ B (P-Pd?dt i<%, 

where pd is the desired pressure. The design problem is 
now treated as a control problem where the control func- 
tion is the wing shape, which is to be chosen to minimize 
I subject to the constraints defined by the flow equations 
(71-72). A variation in the shape will cause a variation 
6p in the pressure and consequently the a variation in the 
cost function 


61 = 



(p-Pd)<5p 


(74) 


Since p depends on w through the equation of state 
(69-70), the variation 6p can be determined from the 
variation 6w. Define the Jacobian matrices 

A, = Ci = JK~- l Aj. (75) 

Then the equation for 6w in the steady state becomes 


— (6Fi) = 0, 

°0 


(76) 


where 


6F{ — CiSw + 6 



fr 


Now, multiplying by a vector co-state variable and 
integrating over the domain 



dij = o, 


and if ip is differentiable this may be integrated by parts 
to give 

X,(i %•*)*>- 

where n t are components of a unit vector normal to the 
boundary. Thus the variation in the cost function may 
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now be written 



On the wing surface Bw , «i = n 3 = 0 and it follows 
from equation (73) that 


CH 


t 


2a: x,y- Plane. 


2b: {.q-Plane. 

Figure 2: Sheared Parabolic Mapping. 



(78) 


Suppose now that is the steady state solution of the 
adjoint equation 


_ r T^_ 

dt • 


0 in D. 


(79) 


At the outer boundary incoming characteristics for ip cor- 
respond to outgoing characteristics for 6w. Consequently, 
as in the two-dimensional case, one can choose boundary 
conditions for ip such that 


riiip T Ci6w — 0. 


Then if the coordinate transformation is such that 8 ( J K ~ 1 ) 
is negligible in the far field, the only remaining boundary 
term is 



Thus by letting ip satisfy the boundary condition, 


(i> 2 


2 , , 2 , d&\ , 

^ f- ip $ "T h ipA-z — J = (p - 

OX 1 OX 2 <7X3 7 


Pd) on Bw, (80) 


we find finally that 


Here x = xi, y = X 2 , z = *3 are the Cartesian coordi- 
nates, and £ and rj+S correspond to parabolic coordinates 
generated by the mapping 


X + iy = x 0 + iyo + ^ a (0 {£ + »(»? + 5)} 2 

at a fixed span station <. x 0 (0 and Vo (0 are the coordi- 
nates of a singular line which is swept to lie just inside the 
leading edge of a swept wing, while a (Q is a scale factor 
to allow for spanwise chord variations. The surface rj = 0 
is a shallow bump corresponding to the wing surface, with 
a height *S(£,0 determined by the equation 

£ + *S= yjl (*b w + *2 /B w )> 


where xb w ( z ) and vb w (*) are coordinates of points ly- 
ing on the wing surface. We now treat S (£, 0 as the 
control. 

In this case the transformation matrix becomes 


K = 


where 


<*(£-(»? + S)$l) -a (rj + S) ,A-a (*7 + 5)S c 
a (*? + & + a£ B + 

0 0 1 

Xjj A + XrjS^ 

Vi Vrj B + VrjSi , 

0 0 1 


A . 

A = a c h*o c , 


B = a c 


r - yo 


+ y<v 



A convenient way to treat a wing is to introduce sheared 
parabolic coordinates as shown in figure 2 through the 
transformation 


* = *o (o + (o iy ~(v+s o) 2 } 

y = yo(0 + a(Oi(T, + S(Z,0) 
z = C 


Now, 

= £ 2 + (»? + 5) 2 
and 


JiC- 1 


y«j Xfj 3 /«jw 4 

-yc *f y?-4 - x f B - J5< 

0 0 J 


Then under a modification 65 

6x( = — a (655^ + (tj + 5)65^) 

6x, = — aSS 

6y ( = a (65 + £6,%) 

6y, = 0. 
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Thus 


6J = 2a 2 (tj + S) 6S 



0 

a6S 

-aB6S ' 

6 ( jr 1 ) = 

-% 

6x£ 

V 


0 

0 

6J 


where 


where A is sufficiently small and non-negative. 

In order to impose a thickness constraint we can define 
a baseline surface So(£,0 below which £(£,<) is not 
allowed to fall. Now if we take A = A (£, 0 as a non- 
negative function such that 

Then the constraint is satisfied, while 


V = 6y{A - Sx^B - a c -6S - <5 JS C - J6S Z . 


Inserting these formulas in equation (81) we find that the 
volume integral in 81 is 


m 


iSf 2 dt dv d( 


- Ill 

+ /// 


tt- {-%/■ + 8x (h + vh) i* d v dc 

dr} 


dtp 7 

dC 


6Jh df] 


where S and 6S are independent of t?. Therefore, inte- 
grating over 7} y the variation in the cost function can be 
reduced to a surface integral of the form 

61 = [ I ( P(M)6S - Q«,C)«S { - R((,06S ( ) d$ dC 

J J B w 

Here 

P — a (tp 2 + S^xp 3 + Ctpi) p 

^ {f/i + (V + S) f 2 + (M + (IJ + S) B) h) dv 


-/ 

-/ 

f w T 

J 


^(fi+Stfi + Ch) dv 


Jd V 

Q = a(^2+(v + S)ip3)p 

+ J {Ui + (» 1 + S)f 2 + (HA + (V + S)B)h} dv 

R = JtpAp 

+ J 

where 


C = 2a(r] + S)S( - A- BS$ + -. 

Also the shape change will be confined to a boundary 
region of the f — C plane, so we can integrate by parts to 
obtain 


61 


f j B vv ^ 




Thus to reduce / we can choose 


ss = -x + 



sg 


&rv 
+ dCj 


d£dt< 0. 


9 IMPLEMENTATION OF THE THREE DIMEN- 
SIONAL METHOD FOR WING DESIGN 

Since three dimensional calculations are much more ex- 
pensive than two dimensional calculations, it is extremely 
important for the practical implementation of the method 
to use fast solution algorithms for the flow and the adjoint 
equations. In this case the author’s FL087 computer pro- 
gram has been used as the basis of the design method. 
FL087 solves the three dimensional Euler equations with 
a cell -centered finite volume scheme, and uses residual 
averaging and multigrid acceleration to obtain very rapid 
steady state solutions, usually in 25 to 50 multigrid cycles 
[8, 9]. Upwind biasing is used to produce nonoscillatory 
solutions, and assure the clean capture of shock waves. 
This is introduced through the addition of carefully con- 
trolled numerical diffusion terms, with a magnitude of 
order Ax 3 in smooth parts of the flow. The program 
corresponds closely to FL082, which was used to imple- 
ment the design method for the two dimensional Euler 
equations. The adjoint equations are treated in the same 
way as the flow equations. The fluxes are first estimated 
by central differences, and then modified by downwind 
biasing through numerical diffusive terms which are sup- 
plied by the same subroutines that were used for the flow 
equations. 

The method has been tested for the optimization of 
a swept wing. The planform was fixed while the wing 
sections were free to be changed arbitrarily by the design 
method. The wing has a unit-semi-span, with 36 degrees 
leading edge sweep. It has a compound trapezoidal plan- 
form, with straight taper from a root chord of 0.38 to a 
chord of 0.26 at the 30 percent span station, and straight 
taper from there to a chord of 0.12 at the tip, with an 
aspect ratio of 8.7. The initial wing sections were based 
on the Kom airfoil, which was designed for shock free 
flow at Mach 0.75 with a lift coefficient of 0.63, and has a 
thickness to chord ratio of 1 1 .5 percent [1]. The thickness 
to chord ratio was increased by a factor of 1 .2 at the root 
and decreased by a ratio of 0.8 at the tip, with a linear vari- 
ation across the span. The inboard sections were rotated 
upwards to give 3.5 degrees twist across the span. 

The two dimensional pressure distribution of the Kom 
airfoil at its design point was introduced as a target pres- 
sure distribution uniformly across the span. This target is 
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presumably not realizable, but serves to favor the estab- 
lishment of relatively benign pressure distribution. The 
total inviscid drag coefficient, due to the combination of 
vortex and shock wave drag, was also included in the cost 
function. Calculations were performed with the lift coef- 
ficient forced to approach a fixed value by adjusting the 
angle of attack every fifth iteration of the flow solution. It 
was found that the computational costs can be reduced by 
using only 15 multigrid cycles in each flow solution, and in 
each adjoint solution. Although this is not enough for full 
convergence, it proves sufficient to provide a shape mod- 
ification which leads to an improvement. Figures 15,16, 
and 17 show the result of a calculation at Mach number of 
0.82, with the lift coefficient forced to approach a value of 
0.5. This calculation was performed on a mesh with 192 
intervals in the £ direction wrapping around the wing, 32 
intervals in the normal rj direction and 48 intervals in the 
spanwise C direction, giving a total of 294912 cells. The 
wing was specified by 33 sections, each with 128 points, 
giving a total of 4224 design variables. The plots show 
the initial wing geometry and pressure distribution, and 
the modified geometry and pressure distribution after 8 
design cycles. The total inviscid drag was reduced from 
0.01 85 to 0.01 1 8. The initial design exhibits a very strong 
shock wave in the inboard region. It can be seen that this 
is completely eliminated, leaving a very weak shock wave 
in the outboard region. The drag reduction is mainly ac- 
complished in the first four design cycles but the pressure 
distribution continues to be adjusted to become more like 
the Kom pressure distribution. 

To verify the solution, the final geometry, after 8 de- 
sign cycles, was analyzed with another method using the 
computer program FL067. This program uses a cell- 
vertex formulation, and has recently been modified to 
incorporate a local extremum diminishing algorithm with 
a very low level of numerical diffusion [12]. When run 
to full convergence it was found that the redesigned wing 
has a drag coefficient of 0.0107 at Mach 0.82 at a lift 
coefficient of 0.5, with a corresponding lift to drag ratio 
of 47. The result is illustrated in Figure 18. A calcu- 
lation at Mach 0.500 shows a drag coefficient of 0.0100 
for a lift coefficient of 0.5. Since in this case the flow is 
entirely subsonic, this provides an estimate of the vortex 
drag for this planform and lift distribution. Thus the de- 
sign method has reduced the shock wave drag coefficient 
to about 0.0007. For a representative transport aircraft the 
parasite drag coefficient of the wing due to skin friction is 
about 0.0050. Also the fuselage drag coefficient is about 
0.0050, the nacelle drag coefficient is about 0.0015, the 
empennage drag coefficient is about 0.0015, and excres- 
cence drag coefficient is about 0.0006. This would give 
a total drag coefficient Cd — 0.0243 for a lift coefficient 
of 0.5, corresponding to a lift to drag ratio L/D = 20.5. 
This would be a substantial improvement over the values 
obtained by currently flying transport aircraft. 

As a further test the redesign was also performed at 
a higher Mach number of 0.85. The initial geometry 


and pressure distributions, and the result of the redesign 
after 10 design cycles are displayed in Figures 19, 20 and 
21. In this case the total inviscid drag was reduced from 
0.026 1 to 0.0 1 32. Again this result has been checked with 
FL067, and when the flow calculation is fully converged, 
it is found that the total inviscid drag coefficient is 0.0118 
at a lift coefficient of 0.5, indicating a shock wave drag 
coefficient of 0.0018. Allowing for the other sources of 
drag for the complete aircraft, it is likely that the best 
operating point for maximum lift to drag ratio would be 
at a somewhat higher lift coefficient. 

10 CONCLUSION 

In the period since this approach to optimal shape design 
was first proposed by the author [ 1 0] , the method has been 
verified by numerical implementation for both potential 
flow and flows modeled by the Euler equations. It has 
been demonstrated that it can be successfully used with a 
finite volume formulation to perform calculations with ar- 
bitrary numerically generated grids [16]. The first results 
which have been obtained for swept wings with the three 
dimensional Euler equations suggest that the method has 
now matured to the point where it can be a very useful 
tool for the design of new airplanes. Even in the case of 
three dimensional flows, the computational requirements 
are so moderate that the calculations can be performed 
with workstations such as the IBM RISC 6000 series. A 
design cycle on a 192x32x48 mesh takes about 1 \ hours 
on an IBM model 530 workstation, allowing overnight 
completion of a design calculation for a swept wing. 
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3a: C p after Zero Design Cycles. 
Design Mach 0.72, C, = 0.5982, C d = 0.0191. 



3b: C p after Zero Design Cycles. 
Design Mach 0.2, C, = 0.5998, C d = -0.0001. 





3c: C p after Eight Design Cycles. 
Design Mach 0.72, C, = 0.5999, C d = 0.0001. 
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3d: C p after Eight Design Cycles. 
Design Mach 0.2, Cj = 0.5998, C d = -0.0001 



Figure 3: Optimization of an Airfoil at Two Design Points. 
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4a: Initial Condition 
C, = 0.0000, C d = 0.0001 




4b: 7 Design Iterations 
Cl = 0.0000, Cd = 0.0000 


Figure 4: Subsonic Non-Lifting Design Case, M = 0.2, a = 0°. 
— , x Initial Airfoil: NACA0012. 

- - -, + Target C p : NACA 64012, M = 0.2. 

Inverse Design 



5a: Initial Condition 
Cl = 0.0000, Cd = 0.0063 


3 

^ _ 





5b: 7 Design Iterations 
C, = 0.0000, Cd = 0.0003 


Figure 5: Transonic Non-Lifting Design Case, M = 0.8 a = 0°. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p : NACA 64021 , M = 0.2. 

Inverse Design 
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6a: Initial Condition 
Cl = 0.0000, Cd = 0.0063 


6b: 8 Design Iterations 
Cl = 0.0000, Cd = 0.0015 


Figure 6: Transonic Non-Lifting Design Case, M = 0.8, a = 0°. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p : NACA 64X, M = 0.8. 

Inverse Design 



7a: Initial Condition 
Ci = 0.7315, C d = 0.0252, a = 2.664* 




7b: 20 Design Iterations 
Ci = 0.7334, C d = 0.0086, a = 0.032° 


Figure 7: Transonic Lifting Design Case, M = 0.735 Fixed Lift. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p \ NACA 64A41 0, M = 0.735, C, = 0.73. 
Inverse Design 
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8a: Initial Condition 8b: 30 Design Iterations 

Cl = 0.5492, Cd = 0.0047, a = 2.709° C, = 0.5496, C d = 0.0045, a = -1.508° 


Figure 8: Transonic Lifting Design Case, M = 0.70, Fixed Lift. 
— , x Initial Airfoil: NACA0012. 

- - -, + Target C p : GAW72, M = 0.70. 

Inverse Design 



9a: Initial Condition 9b: 27 Design Iterations 

Cl = 0.7946, Cd = 0.0358, a = 2.364° Ct = 0.7971, C d = 0.0108, a = 1.053° 


Figure 9: Transonic Lifting Design Case, M = 0.75 Fixed Lift. 
— , x Initial Airfoil: NACA0012. 

- - -, + Target C P : RAE, M = 0.75. 

Inverse Design 
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10a: Initial Condition 
Ct = 0.5037, C d = 0.0127, a = 1.856° 



10b: 2 Design Iterations 
Ci = 0.5042, Cd = 0.0016, a = 1 .990° 


Figure 10: Transonic Lifting Design Case, M = 0.75, Fixed Lift. 
— , x Initial Airfoil: NACA0012. 

Symetric Drag Minimization. 



Figure 11: Transonic Lifting Design Case, M = 0.735 Fixed Lift. 
— , x Initial Airfoil: NACA64A410. 

Camber Only Drag Minimization. 
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12a: Initial Condition 12b: 40 Design Iterations 

Ci = 0.7019, C d = 0.0015, a = 0.266° C; = 0.6612, C d = 0.0136, a = -0.037° 


Figure 12: Lifting Design Case, M = 0.75, Fixed Lift Mode. 
— , x Initial Airfoil: Korn. 

- - -, + Target C p : NACA 64012, M = 0.75. 

Inverse Design 



13a: Initial Condition 13b: 20 Design Iterations 

Ci = 0.6432, C d = 0.0155, a = -0.229° Ct = 0.6297, C d = 0.0077, a = 0.033° 


Figure 13: Lifting Design Case, M = 0.78 Fixed Lift Mode 
Initial Airfoil: Korn. 

Target C p : Korn, M - 0.75. 

Inverse Design 
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14a: Initial Condition 
Ci = 0.6778, C d = 0.0144, a = -0.096° 


* 
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! J * 

14b: 25 Design Iterations 
Ct = 0.6855, Cd = 0.0010, a = -0.722° 


Figure 14: Lifting Design Case, M = 0.75, Fixed Lift Mode. 
Initial Airfoil: NACA64A410. 

Drag Reduction 
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15a: Initial Wing 

Ci = 0.5001, C A - 0.0185, a = -0.958° 


15b: 8 Design Iterations 
Ci = 0.4929, C d = 0.01 18, a = 0. 


Figure 15: Lifting Design Case, M = 0.82, Fixed Lift Mode. 
Drag Reduction 
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UPPER SURFACE PRESSURE 


LOWER SURFACE PRESSURE 


Figure 16: Lifting Design Case, M = 0.82, Fixed Lift Mode. 
Initial Wing: Modfied Korn. 

C L = 0.5001, C D = 0.0185, a = -0.958° 

Drag Reduction 




UPPER SURFACE PRESSURE 


LOWER SURFACE PRESSURE 


Figure 17: Lifting Design Case, M = 0.82, Fixed Lift Mode. 
Design after 8 cycles 
C L = 0.4929, C D = 0.0118, a = 0.172° 

Drag Reduction 
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18c: span station z = 0.50 18d: span station z — 0.75 

Figure 18: FL067 check on redesigned wing. 

M = 0.82, C L = 0.4975, C D = 0.0107, a = 0.200° 
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19a: Initial Wing 

Ci = 0.5033, C d - 0.0261, a = -1.236° 


19b: 10 Design Iterations 
Ci = 0.4956, C d = 0.0132, a = -0.028° 


Figure 19: Lifting Design Case, M = 0.85, Fixed Lift Mode. 
Drag Reduction 
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UPPER SURFACE PRESSURE 


LOWER SURFACE PRESSURE 


Figure 20: Lifting Design Case, M = 0.85, Fixed Lift Mode. 
Initial Airfoil: Modified Korn. 

C L = 0.5033, C D = 0.0261, a = -1.236° 

Drag Reduction 




UPPER SURFACE PRESSURE 


LOWER SURFACE PRESSURE 


Figure 21: Lifting Design Case, M = 0.85, Fixed Lift Mode. 
Design after 10 cycles 
C L = 0.4956, C D = 0.0132, a = -0.028° 

Drag Reduction 
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